Executive Summary

Przeprowadzana analiza w głównej mierze dotyczyła zbioru zawierającego światowe wskaźniki rozwoju zebrane przez organizację Banku Światowego. Dane te obejmują informacje o możliwościach gospodarczych i rozwoju poszczególnych krajów mierzonym przez ponad 100 statystyk. Ponadto do tego zbioru zostały jeszcze dołączone dodatkowe zbiory danych zawierające informacje o kursach wymiany walut, cenach złota, obrocie bitcoinem oraz miesięcznych wynikach S&P Composite.

Pierwszym krokiem po przeczytaniu zbiorów danych było odpowiednie ich przetransformowanie do postaci, która umożliwi sprawne ich połączenie w jeden wspólny zbiór danych. Po tej operacji przystąpiono do czyszczenia danych, usunięcia niektórych cech oraz uzupełnienia wartości pustych.

Następnie skupiono się na szczegółowej analizie wartości atrybutów, gdzie zostały sprawdzone takie cechy jak:

  • Urban population growth (annual %) - cecha szczególna: od 1970 roku do dzisiaj notowany trend spadkowy
  • Survival to age 65 - cecha szczególna: stała różnica utrzymująca się przez wszystkie lata między wartośćią dla kobiet aniżeli dla mężczyzn - na korzyść kobiet
  • CO2 emissions (kt) - cecha szczególna: od 1970 roku do dzisiaj notowany trend wzrostowy

W kolejnym kroku przystąpiono do poszukiwania najbardziej interesujących korelacji w zbiorze danych. Większość ze znalezionych korealcji okazała się być oczywista, przykładowo wzrost populacji silnie skorelowany ze wzrostem populacji kobiet. Jednak poza tymi oczywistymi udało się znaleźć też kilka bardziej interesujących jak np. oczekiwana długość życia obliczona przy urodzenia skorelowana z przeżywalnością 65 roku życia wśród mężczyżn lub światowa emisja metanu skorelowana z liczbą śmierci poniżej 5 roku życia. Dodatkowo zostały jeszcze znalezione cechy, które najbardziej korelują z wartością cen złota dla 4 wybranych krajów - w większości przypadków były to cechy związane z wartością PKB, eksportem oraz importem dóbr i usług.

Ostatni etap pracy obejmował próbę stworzenia regresora, którego zadaniem będzie przewidywanie cen złota. Algorytmem, który zdecydowano się zastosować do tego problemu był Random Forest. Wyniki metryk otrzymane na wcześniej przygotowanym wyniosły 94,69 dolarów dla RMSE oraz 96% dla R^2. W celu poprawy wyników zdecydowano powtórzyć eksperyment na zbiorze stworzonym z pominięciem zbioru “World_Development_Indicators”. Podejście takie przyniosło bardzo obiecujące wyniki na poziomie 8,57 dolara dla RMSE oraz 99% dla R^2.

Analiza ważności atrybutów jednak wykazała, że za tak optymistyczne wyniki odpowiedzialne było przemieszanie przypadków między zbiorem treningowym a testowym (predykcja danych opartych na aspekcie czasowym). Głównym atrybutem mającym największy udział w predykcji był atrybut daty, podczas gdy pozostałe atrybuty charakteryzowały się bardzo małą lub nawet zerową ważnością.

W celu poprawy wyników zdecydowano się na zmianę techniki generowania zbioru treningowego oraz testowego. Metryki dla algorytmu Random Forest dla takiej kombinacji wyniosły 530 dolarów dla RMSE oraz 12% dla R^2. Na nowo utworzonym zbiorze danych zdecydowano się także przeprowadzić analizę z wykorzystaniem algorytmu ARIMA. Wyniki otrzymane dla tego algorytmu są lepsze aniżeli dla Random Forest - wartość RMSE wyniosła około 374 dolary.

Wykorzystane biblioteki

library(tidyverse)
library(readxl)
library(dplyr)
library(ggplot2)
library(kableExtra)
library(zoo)
library(plotly)
library(gganimate)
library(caret)
library(forecast)

Wstępne przetwarzanie zbioru danych

Odczyt danych

W ramach analizy został dostarczony główny zbiór danych stworzony przez Bank Światowy, który zawieraja informacje o możliwościach gospodarczych i rozwoju poszczególnych krajów mierzonym przez ponad 100 statystyk. Oprócz tego, w ramach projektu zostały zebrane informacje o kursach wymiany walut, cenach złota, obrocie bitcoinem oraz miesięcznych wynikach S&P Composite.

indicators_col_names <- c("Country Name", "Country Code", "Series Name", "Series Code", 1970:2020)
indicators_df <- read_xlsx("World_Development_Indicators.xlsx", n_max = 44305, skip = 1, col_names = indicators_col_names)

gold_prices_df <- read.csv("Gold prices.csv")
bc_prices_df <- read.csv("Bitcoin/BCHAIN-MKPRU.csv")
composite_prices_df <- read.csv("S&P Composite.csv")
dim(indicators_df)
## [1] 44304    55

Główny, surowy zbiór danych “World_Development_Indicators” wstępnie zawiera 44304 wiersze oraz 55 kolumn. W kolejnym kroku wspomniany zbiór danych zostanie odpowiednio wyczyszczony w celach dalszej analizy.

Czyszczenie i transformacja wejściowych zbiorów danych

Prezentowana sekcja obejmuje czyszczenie i transformację wejściowych zbiorów danych.

World Development Indicators

Na wstępie spójrzmy na strukturę interesującego nas zbioru danych:

head(indicators_df[, 1:7])
## # A tibble: 6 x 7
##   `Country Name` `Country Code` `Series Name` `Series Code` `1970` `1971` `1972`
##   <chr>          <chr>          <chr>         <chr>         <chr>  <chr>  <chr> 
## 1 Afghanistan    AFG            Urban popula~ SP.URB.GROW   5.748~ 5.860~ 5.899~
## 2 Afghanistan    AFG            Urban popula~ SP.URB.TOTL.~ 11.643 12.021 12.41 
## 3 Afghanistan    AFG            Value lost d~ IC.FRM.OUTG.~ ..     ..     ..    
## 4 Afghanistan    AFG            Urban popula~ SP.URB.TOTL   13009~ 13794~ 14632~
## 5 Afghanistan    AFG            Urban land a~ AG.LND.TOTL.~ ..     ..     ..    
## 6 Afghanistan    AFG            Unemployment~ SL.UEM.TOTL.~ ..     ..     ..

Jak możemy zauważyć mamy tutaj 4 kolumny, które opisują nasze obserwacje (Country Name, Country Code, Series Name, Series Code) oraz pozostałe 51 kolumn, które przechowują wartości odnotowywane dla danej obserwacji w kolejnych latach.

Chcemy oczyścić nasz zbiór danych, dlatego w pierwszym kroku pozbędziemy się kolumn “Country Name” oraz “Series Name”, które nie niosą żadnych dodatkowych informacji, które będą wymagane w dalszej fazie analizy - wartości dla tych kolumn są odpowiednio odwzorowywane przez kolumny “Country Code” oraz “Series Code”

clean_indicators_df <- indicators_df[,c("Country Code", "Series Code", 1970:2020)]

Kolejnym krokiem będzie zamiana wartości “..”, które zostały odczytane z dostarczonego pliku na wartości puste - “NA”. W ten sposób ułatwimy sobie dalszą analizę zbioru

clean_indicators_df[clean_indicators_df == ".."] <- NA

Następnie dokonamy transformacji naszego zbioru danych z wykorzystaniem dplyr’a:

  • Kolumna z etykietami zawierającymi kolejne lata zostanie scalona do jednej nowej kolumny “year”
  • Dla kolumn “Country Code” oraz “Series Code” nastąpi zmiana ich nazw - odpowiednio “c_code” oraz “s_code”
  • Kolumna zawierająca kody wskaźników (“s_code”) zostanie rozszerzona w analizowanym zbiorze danych, tzn. każda wartość wskaźnika stanie się kolumną przechowywującą odpowiednią wartość wskaźnika dla danego kraju i roku

Po wykonaniu tych kroków zbiór danych prezentuje się następująco:

clean_indicators_df <- clean_indicators_df %>%
  pivot_longer(cols = grep(1970, colnames(clean_indicators_df)):grep(2020, colnames(clean_indicators_df)),
               names_to = "year", values_to = "value") %>%
  rename(c_code = "Country Code", s_code = "Series Code") %>%
  pivot_wider(names_from = s_code, values_from = value)

head(clean_indicators_df[, 1:6])
## # A tibble: 6 x 6
##   c_code year  SP.URB.GROW        SP.URB.TOTL.IN.ZS IC.FRM.OUTG.ZS SP.URB.TOTL
##   <chr>  <chr> <chr>              <chr>             <chr>          <chr>      
## 1 AFG    1970  5.7484880426365912 11.643            <NA>           1300949    
## 2 AFG    1971  5.8601019775850567 12.021            <NA>           1379464    
## 3 AFG    1972  5.8992990692598726 12.41             <NA>           1463291    
## 4 AFG    1973  5.8235730871367268 12.809            <NA>           1551037    
## 5 AFG    1974  5.630224006217686  13.219            <NA>           1640869    
## 6 AFG    1975  5.3432279050454463 13.641            <NA>           1730929

Z racji tego, że po transformacji wszystkie kolumny są oznaczone jako typ “character” dokonamy zamiany typów dla wszystkich kolumn poza kolumną “c_code” na typ numeryczny.

clean_indicators_df[, 2:ncol(clean_indicators_df)] <- lapply(clean_indicators_df[, 2:ncol(clean_indicators_df)], as.numeric)

Gold prices

W przypadku zbioru danych zawierającego ceny złota na przestrzeni lat naszym celem będzie transformacja tego zbioru do postaci, dzięki której będziemy mogli połączyć ten zbiór z wcześniej przetwarzanym World Development Indicators. Ostateczny zbiór będzie zawierał następujące kolumny: rok, kurs otwarcia (AM) w dolarach, kurs zamknięcia (PM) w dolarach. Wybrano wydobyć ze zbioru ceny w dolarach w celu zachowania spójności z pozostałymi zbiorami danych. Większość cech w zbiorze World Development Indicators zawiera ceny opisane w dolarach oraz później przetwarzany zbiór danych dotyczący bitcoina także zawiera ceny wyłącznie w dolarach.

Aby przetransformować zbiór do tej postaci należy kolejno:

  • Stworzyć nową kolumnę “year”, dla której wartości wyekstrahujemy z kolumny “Date”
  • Wyfiltrować interesujące nas przedział lat 1970 - 2020
  • Wyfiltrować wiersze, dla których kurs otwarcia i kurs zamknięcia nie jest wartością pustą
  • Pogrupować po nowo stworzonej kolumnie “year”
  • Jako wartości kursu otwarcia i kursu zamknięcia dla danego przyjmujemy pierwsze wartości w grupie - ostatnie w danym roku (zbiór posortowany według dat malejąco)
gold_prices_by_year_df <- gold_prices_df %>%
  mutate(year = format(as.Date(gold_prices_df$Date, format = "%Y-%m-%d"), "%Y")) %>%
  filter(year >= 1970 & year <= 2020) %>%
  filter(!is.na(USD..AM.) & !is.na(USD..PM.)) %>%
  group_by(year) %>%
  summarise(USD_AM = first(USD..AM.), USD_PM = first(USD..PM.))

Dlaczego jako wartości kursu wybieramy pierwsze wartości z grupy, a nie np. średnią czy medianę? Otóż poczynione zostało założenie, że wartości znajdujące się w zbiorze World Development Indicators gromadzone są na stan zakończenia roku kalendarzowego, w związku z czym chcąc zachować spójność między tymi dwoma zbiorami wartości złota zostały wybrane także na stan zakończenia danego roku kalendarzowego.

S&P Composite

W przypadku zbioru zawierającego dane dotyczące S&P Composite transformacja będzie wyglądała analogicznie jak w przypadku zbioru zawierającego ceny złota. Kolumny, które zostały wybrane z tego zbioru to: “S&P Composite” oraz “Real Price”, gdyż uznano je za interesujące oraz mogące mieć wpływ na znalezione zależności we właściwym przetwarzaniu.

composite_prices_by_year_df <- composite_prices_df %>%
  mutate(year = format(as.Date(composite_prices_df$Year, format = "%Y-%m-%d"), "%Y")) %>%
  filter(year >= 1970 & year <= 2020) %>%
  group_by(year) %>%
  summarise(COMPOSITE = first(S.P.Composite), PRICE = first(Real.Price))

Bitcoin

Kolejno przejdziemy do przetworzenia ostatniego zbioru danych wykorzystanego w analizie, który zawiera dane dotyczące Bitcoina. Ze wszystkich dostarczonych zbiorów danych powiązanych z Bitcoinem zdecydowano się na wykorzystanie danych odnośnie ceny rynkowej Bitcoina mierzonej w dolarach (zbiór “BCHAIN-MKPRU”). Wykorzystano jedynie ten zbiór danych z wszystkich dostarczonych zbiorów, gdyż zawiera informacje najbardziej przystępne dla końcowych odbiorców. Transformacja zbioru wygląda analogicznie jak w przypadku zbiorów S&P Composite oraz Gold prices z tą różnicą, że dla zbioru danych dotyczącego Bitcoina posiadamy jedynie informacje dla lat z przedziału 2009 - 2021, dlatego wartości Bitcoina dla lat 1970 - 2008 zostały uzupełnione wartością 0.

missing_bc_prices_df <- data.frame(year = 1970:2008, BC_PRICE = 0.0)
bc_prices_by_year_df <- bc_prices_df %>%
  mutate(year = format(as.Date(bc_prices_df$Date, format = "%Y-%m-%d"), "%Y")) %>%
  filter(year <= 2020) %>%
  group_by(year) %>%
  summarise(BC_PRICE = first(Value)) %>%
  rbind(missing_bc_prices_df)

Przygotowanie i czyszczenie ostatecznego zbioru danych

Prezentowana sekcja obejmuje Przygotowanie i czyszczenie ostatecznego zbioru danych.

Połączenie przygotowanych zbiorów danych

W pierwszym kroku połączymy przygotowane wcześniej zbiory danych w jeden wspólny zbiór danych, który posłuży nam do dalszej analizy.

merged_indicators_df <- clean_indicators_df %>%
  merge(gold_prices_by_year_df, by = "year") %>%
  merge(bc_prices_by_year_df, by = "year") %>%
  merge(composite_prices_by_year_df, by = "year")

Czyszczenie połączonego zbioru danych

Przygotowany w poprzednim punkcie zbiór danych nadal zawiera wiele wartości pustych, którymi zajmiemy się w następnej kolejności.

sort(colSums(is.na(merged_indicators_df)), decreasing = T) %>%
  kable %>%
  kable_styling("striped") %>%
  scroll_box(height = "250px")
x
IC.FRM.OUTG.ZS 10314
IC.TAX.METG 10312
SH.STA.DIAB.ZS 10207
DT.ODA.OATL.CD 10190
FX.OWN.TOTL.ZS 10174
AG.LND.TOTL.UR.K2 10086
SE.SEC.TCAQ.UP.ZS 10044
EN.POP.SLUM.UR.ZS 9983
SE.SEC.TCAQ.ZS 9818
NY.GSR.NFCY.KN 9742
SH.ALC.PCAP.LI 9714
IC.BUS.DFRN.XQ 9686
SE.ADT.LITR.ZS 9464
SM.POP.TOTL.ZS 9386
IC.LGL.CRED.XQ 9318
SE.PRM.TCAQ.ZS 9294
SL.UEM.NEET.MA.ZS 9185
SL.UEM.NEET.FE.ZS 9184
SL.UEM.NEET.ZS 9165
SP.POP.SCIE.RD.P6 9062
AG.LND.PRCP.MM 9016
FB.BNK.CAPA.ZS 9007
ER.H2O.INTR.PC 8944
ER.H2O.INTR.K3 8939
SI.DST.10TH.10 8909
SL.TLF.PART.ZS 8905
SL.UEM.ADVN.ZS 8849
DT.NFL.BOND.CD 8820
SE.SEC.ENRL.UP.TC.ZS 8720
GB.XPD.RSDV.GD.ZS 8675
IS.RRS.TOTL.KM 8662
IC.ELC.TIME 8657
IS.RRS.PASG.KM 8653
CM.MKT.INDX.ZG 8646
IS.RRS.GOOD.MT.K6 8552
IT.NET.SECR 8444
EN.ATM.PM25.MC.T1.ZS 8424
EN.ATM.PM25.MC.T2.ZS 8424
EN.ATM.PM25.MC.T3.ZS 8424
IT.NET.SECR.P6 8370
GC.AST.TOTL.GD.ZS 8357
EN.ATM.PM25.MC.M3 8340
EN.ATM.PM25.MC.ZS 8340
CM.MKT.TRNR 8201
IC.WRH.DURS 8028
IC.TAX.PAYM 7990
CM.MKT.TRAD.GD.ZS 7882
FB.ATM.TOTL.P5 7868
CM.MKT.TRAD.CD 7849
IC.LGL.DURS 7713
SN.ITK.DEFC.ZS 7606
TX.VAL.ICTG.ZS.UN 7594
GC.TAX.GSRV.VA.ZS 7433
SH.XPD.CHEX.GD.ZS 7311
SH.XPD.CHEX.PC.CD 7178
SH.STA.SUIC.P5 7028
SH.STA.SUIC.FE.P5 7028
SH.STA.SUIC.MA.P5 7028
SH.STA.TRAF.P5 7020
IP.JRN.ARTC.SC 6998
SE.TER.ENRL.TC.ZS 6956
GC.XPN.TOTL.GD.ZS 6939
ST.INT.XPND.CD 6889
GC.TAX.EXPT.ZS 6888
GC.XPN.INTP.ZS 6882
GC.TAX.EXPT.CN 6875
IP.PAT.RESD 6868
GC.TAX.GSRV.CN 6800
SE.XPD.TOTL.GD.ZS 6787
GC.TAX.YPKG.ZS 6774
GC.TAX.YPKG.CN 6768
GC.TAX.INTT.CN 6764
EN.ATM.GHGT.ZG 6762
GC.TAX.TOTL.CN 6737
GC.TAX.GSRV.RV.ZS 6723
FR.INR.RINR 6703
GC.TAX.YPKG.RV.ZS 6683
GC.TAX.INTT.RV.ZS 6683
GC.TAX.TOTL.GD.ZS 6676
DT.DOD.DSTC.XP.ZS 6640
IP.PAT.NRES 6601
FR.INR.LEND 6550
SG.GEN.PARL.ZS 6404
FR.INR.DPST 6397
IP.TMK.RESD 6357
DT.DOD.DSTC.IR.ZS 6339
IP.TMK.NRES 6326
EN.ATM.METH.ZG 6296
EN.ATM.NOXE.ZG 6274
SE.PRE.ENRL.TC.ZS 6197
SL.UEM.TOTL.NE.ZS 6194
SE.SEC.ENRL.TC.ZS 6160
IP.TMK.TOTL 6083
DT.DOD.DECT.GN.ZS 5961
SE.SEC.TCHR 5894
SE.ENR.TERT.FM.ZS 5777
EN.ATM.CO2E.PP.GD.KD 5714
EN.ATM.CO2E.PP.GD 5615
DT.DOD.DSTC.ZS 5508
NY.GNP.MKTP.KD.ZG 5393
SL.EMP.SELF.MA.ZS 5301
SL.EMP.SELF.ZS 5301
SL.EMP.SELF.FE.ZS 5301
SL.IND.EMPL.ZS 5301
SL.SRV.EMPL.ZS 5301
SL.AGR.EMPL.ZS 5301
SL.EMP.MPYR.ZS 5301
EG.ELC.ACCS.ZS 5253
EG.ELC.RNEW.ZS 5196
NY.GNS.ICTR.CD 5122
SL.TLF.TOTL.IN 4938
NY.GNS.ICTR.ZS 4902
NY.TAX.NIND.KN 4856
EN.ATM.CO2E.EG.ZS 4855
EN.CO2.TRAN.ZS 4839
EN.CO2.BLDG.ZS 4839
EN.CO2.OTHX.ZS 4839
EN.CO2.MANF.ZS 4839
EN.CO2.ETOT.ZS 4839
EG.ELC.NUCL.ZS 4835
IT.NET.USER.ZS 4786
EG.FEC.RNEW.ZS 4776
EG.ELC.RNWX.ZS 4731
EG.ELC.FOSL.ZS 4731
EG.ELC.COAL.ZS 4731
EG.ELC.HYRO.ZS 4731
EG.ELC.NGAS.ZS 4731
EG.ELC.RNWX.KH 4725
TX.VAL.FUEL.ZS.UN 4666
SE.PRM.ENRL.TC.ZS 4628
EN.URB.MCTY 4590
BX.PEF.TOTL.CD.WD 4579
BN.KLT.PTXL.CD 4436
NE.EXP.GNFS.KD.ZG 4428
TX.VAL.FOOD.ZS.UN 4301
TM.VAL.FUEL.ZS.UN 4275
TM.VAL.FOOD.ZS.UN 4267
NY.TAX.NIND.CD 4116
BN.GSR.FCTY.CD 4105
TX.VAL.TRAN.ZS.WT 4094
NY.TAX.NIND.CN 4072
BG.GSR.NFSV.GD.ZS 4044
NV.IND.MANF.ZS 3956
TM.VAL.TRAN.ZS.WT 3931
BX.GSR.FCTY.CD 3912
BM.GSR.FCTY.CD 3905
SE.SEC.ENRL 3890
BX.GSR.MRCH.CD 3850
BM.GSR.NFSV.CD 3841
BM.GSR.MRCH.CD 3840
BX.GSR.NFSV.CD 3833
NV.SRV.TOTL.ZS 3821
DT.ODA.ODAT.CD 3739
NE.DAB.TOTL.ZS 3714
FM.AST.DOMS.CN 3622
FP.CPI.TOTL 3577
FP.CPI.TOTL.ZG 3398
NE.DAB.TOTL.CD 3396
NY.GDS.TOTL.CD 3351
NY.GDS.TOTL.ZS 3305
EN.ATM.NOXE.EG.ZS 3270
EN.ATM.CO2E.KD.GD 3095
EN.URB.LCTY 3060
NY.GSR.NFCY.CD 2991
NE.TRD.GNFS.ZS 2957
NE.IMP.GNFS.ZS 2956
NY.GSR.NFCY.CN 2952
NE.EXP.GNFS.CD 2950
NE.IMP.GNFS.CD 2949
EN.URB.LCTY.UR.ZS 2726
NY.GDP.NGAS.RT.ZS 2653
EN.ATM.CO2E.SF.ZS 2485
EN.ATM.CO2E.LF.ZS 2485
EN.ATM.CO2E.GF.ZS 2485
TX.VAL.MRCH.HI.ZS 2228
EN.ATM.CO2E.SF.KT 2159
EN.ATM.CO2E.GF.KT 2159
NY.GDP.PCAP.KD.ZG 2114
EN.ATM.CO2E.LF.KT 2112
NY.GDP.MKTP.KD.ZG 2111
SH.DTH.MORT 2094
NY.GDP.TOTL.RT.ZS 2009
EN.ATM.CO2E.PC 1997
EN.ATM.CO2E.KT 1994
NY.GDP.PCAP.CD 1859
NY.GDP.MKTP.CD 1856
SP.DYN.IMRT.IN 1818
EN.ATM.GHGT.KT.CE 1556
EN.ATM.METH.KT.CE 1374
EN.ATM.NOXE.KT.CE 1354
EN.ATM.METH.EG.KT.CE 1112
SP.DYN.TO65.FE.ZS 1108
SP.DYN.TO65.MA.ZS 1108
SP.DYN.LE00.IN 1022
SP.POP.TOTL.MA.IN 950
SP.POP.TOTL.FE.IN 950
SP.POP.TOTL.MA.ZS 927
SP.POP.TOTL.FE.ZS 927
SP.POP.65UP.TO.ZS 927
SP.POP.1564.TO.ZS 927
SP.POP.0014.TO.ZS 927
SE.PRM.AGES 833
SP.DYN.CBRT.IN 815
SP.RUR.TOTL.ZG 462
EN.POP.DNST 139
AG.LND.TOTL.K2 116
SP.URB.GROW 85
SP.URB.TOTL 83
SP.RUR.TOTL 83
SP.URB.TOTL.IN.ZS 60
SP.RUR.TOTL.ZS 60
SP.POP.GROW 35
SP.POP.TOTL 32
year 0
c_code 0
USD_AM 0
USD_PM 0
BC_PRICE 0
COMPOSITE 0
PRICE 0

 

W pierwszym kroku usuniemy kolumny, które zawierają więcej niż 30% wartości pustych. Zdecydowano się na taki krok, ponieważ uznano, że takie kolumny nie przyniosą wartości dodanej dalszej analizie - nie zostaną znalezione interesujące korelacje z wykorzystaniem takich kolumn oraz takie cechy nie będą miały dużego wkładu przy problemie uczenia maszynowego ze względu na istniejące braki.

merged_indicators_df <- merged_indicators_df[, which(colMeans(is.na(merged_indicators_df)) <= 0.3)]

Dla kolumn, które pozostały w zbiorze po filtrowaniu przeprowadzimy operację uzupełnienia wartości pustych zgodnie z przyjętą strategią:

  • Wartości puste zostaną uzupełnione średnią dla danego kraju biorąc jedynie pod uwagę rok poprzedni oraz rok następny
  • Wartości puste na krańcach dla danej cechy i dla danego kraju zostaną uzupełnione najbliższymi wartościami niepustymi, przykładowo dla wektora [NA, 2, 5, 7, NA] jego postać po transformacji będzie wyglądała następująco [2, 2, 5, 7, 7]
clean_df <- merged_indicators_df %>%
  group_by(c_code) %>%
  mutate_at(vars(-group_cols()), na.approx, na.rm = FALSE, rule = 2) %>%
  ungroup()

Po wykonaniu tej transformacji otrzymujemy oczyszczony, ostateczny zbiór danych, który wykorzystamy w dalszej analizie.

Podsumowanie przygotowanego zbioru danych oraz podstawowe statystyki

year c_code SP.URB.GROW SP.URB.TOTL.IN.ZS SP.URB.TOTL NE.TRD.GNFS.ZS NY.GDP.TOTL.RT.ZS EN.ATM.GHGT.KT.CE SP.DYN.TO65.FE.ZS SP.DYN.TO65.MA.ZS SP.RUR.TOTL.ZG SP.RUR.TOTL.ZS SP.RUR.TOTL SE.PRM.AGES SP.POP.TOTL SP.POP.TOTL.MA.IN SP.POP.TOTL.MA.ZS SP.POP.TOTL.FE.ZS SP.POP.TOTL.FE.IN EN.URB.LCTY.UR.ZS EN.URB.LCTY SP.POP.GROW EN.POP.DNST SP.POP.65UP.TO.ZS SP.POP.1564.TO.ZS SP.POP.0014.TO.ZS SH.DTH.MORT EN.ATM.NOXE.KT.CE NY.GSR.NFCY.CD NY.GSR.NFCY.CN NY.GDP.NGAS.RT.ZS SP.DYN.IMRT.IN EN.ATM.METH.KT.CE EN.ATM.METH.EG.KT.CE TX.VAL.MRCH.HI.ZS SP.DYN.LE00.IN AG.LND.TOTL.K2 NE.IMP.GNFS.CD NE.IMP.GNFS.ZS NY.GDP.PCAP.CD NY.GDP.PCAP.KD.ZG NY.GDP.MKTP.KD.ZG NY.GDP.MKTP.CD NE.EXP.GNFS.CD EN.ATM.CO2E.SF.ZS EN.ATM.CO2E.SF.KT EN.ATM.CO2E.LF.KT EN.ATM.CO2E.LF.ZS EN.ATM.CO2E.GF.KT EN.ATM.CO2E.GF.ZS EN.ATM.CO2E.PC EN.ATM.CO2E.KT EN.ATM.CO2E.KD.GD SP.DYN.CBRT.IN USD_AM USD_PM BC_PRICE COMPOSITE PRICE
Min. :1970 Length:10608 Min. :-187.142 Min. : 2.845 Min. :1.267e+03 Min. : 0.021 Min. : 0.0000 Min. : 1 Min. : 6.464 Min. : 1.477 Min. :-235.7924 Min. : 0.00 Min. :0.000e+00 Min. :4.000 Min. :5.740e+03 Min. :2.528e+04 Min. :44.37 Min. :23.29 Min. :2.586e+04 Min. : 2.867 Min. : 18587 Min. :-10.9551 Min. : 0.136 Min. : 0.6856 Min. :45.45 Min. :11.05 Min. : 0 Min. : 0.0 Min. :-9.905e+10 Min. :-4.813e+14 Min. : 0.0000 Min. : 1.50 Min. : 0 Min. : 0 Min. : 0.0074 Min. :18.91 Min. : 2 Min. :0.000e+00 Min. : 0.00 Min. : 22.8 Min. :-64.9924 Min. :-64.047 Min. :8.824e+06 Min. :6.933e+05 Min. : -4.324 Min. : -114 Min. : -161 Min. : -6.089 Min. : -147 Min. : -0.7295 Min. : 0.0000 Min. : 0 Min. :0.0000 Min. : 5.90 Min. : 37.38 Min. : 37.38 Min. : 0 Min. : 67.07 Min. : 352.8
1st Qu.:1982 Class :character 1st Qu.: 1.036 1st Qu.: 33.380 1st Qu.:3.464e+05 1st Qu.: 47.375 1st Qu.: 0.1327 1st Qu.: 5228 1st Qu.:60.747 1st Qu.:52.028 1st Qu.: -0.4820 1st Qu.:26.15 1st Qu.:2.781e+05 1st Qu.:6.000 1st Qu.:7.843e+05 1st Qu.:9.653e+05 1st Qu.:48.96 1st Qu.:49.63 1st Qu.:9.448e+05 1st Qu.: 20.465 1st Qu.: 610060 1st Qu.: 0.6114 1st Qu.: 23.407 1st Qu.: 3.2328 1st Qu.:53.43 1st Qu.:23.38 1st Qu.: 681 1st Qu.: 378.6 1st Qu.:-7.350e+08 1st Qu.:-7.668e+09 1st Qu.: 0.0000 1st Qu.: 12.90 1st Qu.: 1320 1st Qu.: 50 1st Qu.: 55.8564 1st Qu.:59.92 1st Qu.: 16262 1st Qu.:1.046e+09 1st Qu.: 25.90 1st Qu.: 712.2 1st Qu.: -0.4164 1st Qu.: 1.075 1st Qu.:1.982e+09 1st Qu.:6.980e+08 1st Qu.: 0.000 1st Qu.: 0 1st Qu.: 803 1st Qu.: 40.242 1st Qu.: 0 1st Qu.: 0.0000 1st Qu.: 0.5582 1st Qu.: 1041 1st Qu.:0.2403 1st Qu.:14.64 1st Qu.: 290.85 1st Qu.: 290.20 1st Qu.: 0 1st Qu.: 139.40 1st Qu.: 560.1
Median :1995 Mode :character Median : 2.365 Median : 53.132 Median :2.361e+06 Median : 70.669 Median : 1.6701 Median : 29270 Median :77.192 Median :64.691 Median : 0.6466 Median :46.87 Median :2.244e+06 Median :6.000 Median :5.218e+06 Median :3.217e+06 Median :49.66 Median :50.34 Median :3.245e+06 Median : 30.890 Median : 1287166 Median : 1.5462 Median : 69.545 Median : 4.7232 Median :60.61 Median :34.25 Median : 5550 Median : 3115.0 Median :-7.079e+07 Median :-2.317e+08 Median : 0.0000 Median : 31.60 Median : 7084 Median : 798 Median : 72.9729 Median :69.27 Median : 107160 Median :4.191e+09 Median : 37.37 Median : 2328.2 Median : 1.9730 Median : 3.609 Median :9.046e+09 Median :3.431e+09 Median : 2.651 Median : 150 Median : 4177 Median : 68.498 Median : 33 Median : 0.7463 Median : 2.4278 Median : 9153 Median :0.3874 Median :23.36 Median : 391.50 Median : 391.75 Median : 0 Median : 614.57 Median :1093.0
Mean :1995 NA Mean : 2.650 Mean : 53.789 Mean :5.560e+07 Mean : 84.643 Mean : 6.4254 Mean : 724068 Mean :71.987 Mean :62.460 Mean : -0.8912 Mean :46.21 Mean :6.936e+07 Mean :6.151 Mean :1.244e+08 Mean :6.865e+07 Mean :49.92 Mean :50.08 Mean :6.748e+07 Mean : 33.786 Mean : 2947763 Mean : 1.6705 Mean : 352.031 Mean : 6.8287 Mean :60.07 Mean :33.10 Mean : 259903 Mean : 55329.9 Mean :-4.078e+08 Mean :-1.486e+12 Mean : 0.2571 Mean : 45.39 Mean : 148724 Mean : 48883 Mean : 67.9698 Mean :66.52 Mean : 2692296 Mean :1.645e+11 Mean : 45.37 Mean : 9834.2 Mean : 1.6456 Mean : 3.334 Mean :6.135e+11 Mean :1.679e+11 Mean : 16.491 Mean : 211197 Mean : 170212 Mean : 65.679 Mean : 88086 Mean : 12.2118 Mean : 4.9492 Mean : 493766 Mean :0.5985 Mean :25.98 Mean : 590.16 Mean : 590.32 Mean : 1108 Mean : 901.85 Mean :1299.5
3rd Qu.:2008 NA 3rd Qu.: 3.956 3rd Qu.: 73.850 3rd Qu.:9.242e+06 3rd Qu.:101.847 3rd Qu.: 7.7965 3rd Qu.: 107835 3rd Qu.:84.711 3rd Qu.:74.050 3rd Qu.: 1.8399 3rd Qu.:66.62 3rd Qu.:9.390e+06 3rd Qu.:7.000 3rd Qu.:1.936e+07 3rd Qu.:1.147e+07 3rd Qu.:50.37 3rd Qu.:51.04 3rd Qu.:1.155e+07 3rd Qu.: 43.017 3rd Qu.: 2996023 3rd Qu.: 2.5844 3rd Qu.: 164.397 3rd Qu.: 9.8936 3rd Qu.:66.10 3rd Qu.:43.20 3rd Qu.: 42858 3rd Qu.: 11881.7 3rd Qu.: 1.620e+06 3rd Qu.: 3.200e+06 3rd Qu.: 0.0797 3rd Qu.: 67.78 3rd Qu.: 28394 3rd Qu.: 5541 3rd Qu.: 85.0712 3rd Qu.:74.27 3rd Qu.: 553185 3rd Qu.:2.357e+10 3rd Qu.: 56.27 3rd Qu.: 10407.9 3rd Qu.: 4.3091 3rd Qu.: 5.924 3rd Qu.:6.409e+10 3rd Qu.:2.395e+10 3rd Qu.: 26.986 3rd Qu.: 8911 3rd Qu.: 29378 3rd Qu.: 94.248 3rd Qu.: 9604 3rd Qu.: 18.7212 3rd Qu.: 6.8254 3rd Qu.: 61910 3rd Qu.:0.7018 3rd Qu.:36.68 3rd Qu.: 870.00 3rd Qu.: 869.75 3rd Qu.: 0 3rd Qu.:1330.93 3rd Qu.:1916.2
Max. :2020 NA Max. : 48.936 Max. :100.000 Max. :4.352e+09 Max. :860.800 Max. :87.5075 Max. :45873850 Max. :96.093 Max. :92.978 Max. : 29.6283 Max. :97.16 Max. :3.399e+09 Max. :8.000 Max. :7.753e+09 Max. :3.907e+09 Max. :76.71 Max. :55.63 Max. :3.843e+09 Max. :100.000 Max. :37468302 Max. : 17.6334 Max. :21388.600 Max. :28.3973 Max. :86.40 Max. :51.57 Max. :12493789 Max. :2986520.0 Max. : 2.923e+11 Max. : 1.051e+14 Max. :22.4135 Max. :219.30 Max. :8174420 Max. :3187680 Max. :100.0000 Max. :85.42 Max. :129956634 Max. :2.472e+13 Max. :427.58 Max. :190512.7 Max. :140.3670 Max. :149.973 Max. :8.761e+13 Max. :2.525e+13 Max. :216.648 Max. :15291329 Max. :10482498 Max. :258.524 Max. :7056781 Max. :207.3675 Max. :360.8532 Max. :34041046 Max. :5.3510 Max. :56.95 Max. :1877.55 Max. :1887.60 Max. :28857 Max. :3695.31 Max. :3873.1
NA NA NA’s :51 NA’s :51 NA’s :51 NA’s :765 NA’s :153 NA’s :510 NA’s :918 NA’s :918 NA’s :357 NA’s :51 NA’s :51 NA’s :612 NA NA’s :918 NA’s :918 NA’s :918 NA’s :918 NA’s :2703 NA’s :3060 NA NA NA’s :918 NA’s :918 NA’s :918 NA’s :1122 NA’s :408 NA’s :867 NA’s :867 NA’s :612 NA’s :1122 NA’s :408 NA NA’s :510 NA’s :457 NA NA’s :816 NA’s :765 NA’s :153 NA’s :204 NA’s :204 NA’s :153 NA’s :816 NA’s :663 NA’s :459 NA’s :408 NA’s :663 NA’s :459 NA’s :663 NA’s :663 NA’s :663 NA’s :1171 NA’s :204 NA NA NA NA NA

 

Jak możemy zauważyć mimo uzupełnienia wartości pustych to w naszym zbiorze nadal występują cechy, które takie wartości posiadają. Spowodowane jest to tym, że wartości NA dla danej cechy były uzupełniane w obrębie określonego kraju - w sytuacji, gdy dla danego kraju oraz dla danej cechy nie odnotowaliśmy żadnej obserwacji to algorytm uzupełniający wartości puste nie miał na podstawie czego obliczyć pozostałych wartości. Jak możemy zauważyć w wierszu zawierającym wartości puste występują jedynie wielokrotności liczby 51 - wartość ta jest równa liczbie obserwacji dla danego kraju dla danej cechy w latach 1970 - 2020.

dim(clean_df)
## [1] 10608    59

Ostatecznie po czyszczeniu i wszystkich transformacjach nasz zbiór danych zawiera 10608 wierszy oraz 59 kolumn.

Szczegółowa analiza zbioru danych

Prezentowana sekcja obejmuje szczegółową analizę stworzonego w poprzednim punkcie zbioru danych.

Analiza wartości atrybutów

Z racji, że analizowany zbiór danych zawiera 57 atrybutów opisujących obserwację i analiza każdego atrybutu z osobna zajęłaby znaczną część tego raportu, dlatego zdecydowano się przeanalizować kilka interesujących, wybranych cech:

  • Urban population growth (annual %) - SP.URB.GROW
  • Survival to age 65, female (% of cohort) - SP.DYN.TO65.FE.ZS
  • Survival to age 65, male (% of cohort) - SP.DYN.TO65.MA.ZS
  • CO2 emissions (kt) - EN.ATM.CO2E.KT

Urban population growth (annual %)

year SP.URB.GROW
Min. :1970 Min. :-187.142
1st Qu.:1982 1st Qu.: 1.036
Median :1995 Median : 2.365
Mean :1995 Mean : 2.650
3rd Qu.:2008 3rd Qu.: 3.956
Max. :2020 Max. : 48.936
NA NA’s :51

 

growth_df %>%
  group_by(year) %>%
  summarise(growth = mean(SP.URB.GROW, na.rm = T)) %>% 
  ggplot(aes(year, growth, color = "red")) + geom_point() + geom_line() + theme(legend.position = "none")

growth_df <- growth_df %>%
  filter(SP.URB.GROW > -10 & SP.URB.GROW < 20)

interval_col <- findInterval(growth_df$year, seq(min(growth_df$year), max(growth_df$year), 10), rightmost.closed = T)

growth_plot <- growth_df %>%
  cbind(interval_col) %>%
  ggplot(aes(interval_col, SP.URB.GROW, group = interval_col)) + geom_boxplot()

ggplotly(growth_plot)

Wartości dla każdego roku zostały uśrednione, żeby móc odpowiednio przedstawić je za pomocą wizualizacji.

Jak możemy zauważyć na wykresach wartość wzrostu liczby ludności mieszkającej w miastach (mierzonej w %) notuje tendencję malejącą. Najwyższa wartość wystąpiła w roku 1971 i wyniosła 4,13%, najniższa wartość została odnotowana dla roku 2020 i ostatecznie wyniosła 1,78%. Warto także zwrócić uwagę na przypadek odnotowany w 1975, gdzie wystąpił znaczący spadek - główny wpływ na tę wartość miała tocząca się w tamtych czasach wojna domowa w Kambodży.

Na wykresie pudełkowym obserwacje zostały pogrupowane co 10 lat (1970-1979, 1980-1989, …, 2010-2020). Jak możemy zwrócić uwagę środkowa linia pudełka - mediana stopniowo spadała rozpoczynając od wartości 3,46 w pierwszej grupie i kończać na wartości 1,71 w ostatniej grupie. Warto także zwrócić uwagę na znaczą liczbę outlierów dla przedostatniej grupy, co świadczy o zanotowanej tendencji wzrostowej w niektórych krajach.

Survival to age 65

year SP.DYN.TO65.FE.ZS SP.DYN.TO65.MA.ZS
Min. :1970 Min. : 6.464 Min. : 1.477
1st Qu.:1982 1st Qu.:60.747 1st Qu.:52.028
Median :1995 Median :77.192 Median :64.691
Mean :1995 Mean :71.987 Mean :62.460
3rd Qu.:2008 3rd Qu.:84.711 3rd Qu.:74.050
Max. :2020 Max. :96.093 Max. :92.978
NA NA’s :918 NA’s :918

 

survival_df %>%
  group_by(year) %>%
  summarise(female = mean(SP.DYN.TO65.FE.ZS, na.rm = T), male = mean(SP.DYN.TO65.MA.ZS, na.rm = T)) %>%
  pivot_longer(cols = c("female", "male"), names_to = "gender", values_to = "value (%)") %>%
  ggplot(aes(year, `value (%)`, color = gender)) + geom_point() + geom_line()

Wartości przeżywalności 65 roku życia kobiet i mężczyzn notują tend wzrostowy - w miarę upływu lat przeżywalność wśród kobiet i mężczyzn stopniowo wzrasta. Można także zwrócić uwagę na fakt, że przeżywalność wśród kobiet jest znacząco wyższa niż u mężczyzn, co potwierdza znany fakt, że to kobiety żyją dłużej od mężczyzn. Najniższy pojedynczy, odnotowany % przeżywalności u kobiet wyniósł 6,46% natomiast u mężczyzn było to 1,47%. Mediana w przypadku kobiet wyniosła 71,987%, u mężczyzn 64,691%. Najwyższy pojedynczy odnotowany % u kobiet wyniósł 96,093%, u mężczyzn 92,978%.

CO2 emissions (kt)

year EN.ATM.CO2E.KT
Min. :1970 Min. : 0
1st Qu.:1982 1st Qu.: 1041
Median :1995 Median : 9153
Mean :1995 Mean : 493766
3rd Qu.:2008 3rd Qu.: 61910
Max. :2020 Max. :34041046
NA NA’s :663

 

emission_df %>%
  group_by(year) %>%
  summarise(emission = mean(EN.ATM.CO2E.KT, na.rm = T)) %>% 
  ggplot(aes(year, emission, color = "red")) + geom_point() + geom_line() + theme(legend.position = "none")

Wartośc emisji CO2 od roku 1970 do roku 2020 notuje tend wzrostowy. Warty odnotowania jest także fakt znaczącego spadku emisji między rokiem 1989 a rokiem 1990. Mediana odnotowana dla tej cechy wyniosła 9153 kiloton, natomiast najwyższa wartość, która została odnotowana wyniosła 34041046 kiloton.

Badanie korelacji między zmiennymi

W omawianej sekcji przedstawione zostaną wartości współczynnika korelacji Pearsona znalezione w analizowanym zbiorze danych, dla których wartość bezwzględna współczynnika jest większa jak 0,75.

rowname variable correlation
SP.RUR.TOTL.ZS SP.URB.TOTL.IN.ZS -1.0000000
SP.POP.TOTL.FE.ZS SP.POP.TOTL.MA.ZS -1.0000000
SP.POP.TOTL.FE.IN SP.POP.TOTL 0.9999817
SP.POP.TOTL.MA.IN SP.POP.TOTL 0.9999812
SP.POP.TOTL.FE.IN SP.POP.TOTL.MA.IN 0.9999259
NE.EXP.GNFS.CD NE.IMP.GNFS.CD 0.9993788
EN.ATM.CO2E.KT EN.ATM.GHGT.KT.CE 0.9923771
EN.ATM.METH.KT.CE EN.ATM.NOXE.KT.CE 0.9903804
SP.DYN.LE00.IN SP.DYN.TO65.FE.ZS 0.9880276
EN.ATM.METH.KT.CE SP.POP.TOTL.FE.IN 0.9878795
SP.POP.TOTL.MA.IN SP.RUR.TOTL 0.9872349
NY.GDP.MKTP.CD NE.IMP.GNFS.CD 0.9871843
EN.ATM.METH.KT.CE SP.POP.TOTL 0.9871280
SP.POP.TOTL SP.RUR.TOTL 0.9868278
SP.POP.TOTL.FE.IN SP.RUR.TOTL 0.9862969
EN.ATM.METH.KT.CE SP.POP.TOTL.MA.IN 0.9862476
NE.EXP.GNFS.CD NY.GDP.MKTP.CD 0.9840261
EN.ATM.METH.EG.KT.CE EN.ATM.METH.KT.CE 0.9824191
EN.ATM.NOXE.KT.CE EN.ATM.GHGT.KT.CE 0.9820760
AG.LND.TOTL.K2 EN.ATM.NOXE.KT.CE 0.9805704
SP.POP.TOTL.FE.IN SP.URB.TOTL 0.9801023
AG.LND.TOTL.K2 EN.ATM.METH.KT.CE 0.9800440
EN.ATM.METH.EG.KT.CE SP.URB.TOTL 0.9799983
SP.POP.TOTL SP.URB.TOTL 0.9796048
EN.ATM.CO2E.KT EN.ATM.CO2E.SF.KT 0.9795565
SP.POP.TOTL.MA.IN SP.URB.TOTL 0.9789376
EN.ATM.METH.KT.CE SP.URB.TOTL 0.9747790
EN.ATM.CO2E.SF.KT EN.ATM.GHGT.KT.CE 0.9733308
EN.ATM.NOXE.KT.CE SP.POP.TOTL.FE.IN 0.9725662
EN.ATM.CO2E.KT EN.ATM.CO2E.LF.KT 0.9722282
EN.ATM.CO2E.KT EN.ATM.CO2E.GF.KT 0.9721656
EN.ATM.METH.KT.CE EN.ATM.GHGT.KT.CE 0.9721062
EN.ATM.NOXE.KT.CE SP.POP.TOTL 0.9715543
EN.ATM.NOXE.KT.CE SP.POP.TOTL.MA.IN 0.9703351
EN.ATM.METH.EG.KT.CE EN.ATM.GHGT.KT.CE 0.9690934
EN.ATM.CO2E.SF.KT EN.ATM.METH.EG.KT.CE 0.9679511
EN.ATM.METH.KT.CE SP.RUR.TOTL 0.9677979
EN.ATM.CO2E.LF.KT EN.ATM.GHGT.KT.CE 0.9664840
EN.ATM.GHGT.KT.CE SP.URB.TOTL 0.9663866
EN.ATM.METH.EG.KT.CE SP.POP.TOTL.FE.IN 0.9661306
EN.ATM.CO2E.GF.KT EN.ATM.CO2E.LF.KT 0.9653635
EN.ATM.METH.EG.KT.CE SP.POP.TOTL 0.9652718
SP.DYN.LE00.IN SP.DYN.TO65.MA.ZS 0.9651897
EN.ATM.METH.EG.KT.CE SP.POP.TOTL.MA.IN 0.9641086
EN.ATM.CO2E.SF.KT SP.URB.TOTL 0.9640689
EN.ATM.METH.EG.KT.CE EN.ATM.NOXE.KT.CE 0.9632798
EN.ATM.NOXE.KT.CE SP.URB.TOTL 0.9615191
EN.ATM.CO2E.GF.KT EN.ATM.GHGT.KT.CE 0.9575266
EN.ATM.CO2E.KT SP.URB.TOTL 0.9566841
EN.ATM.CO2E.KT EN.ATM.METH.EG.KT.CE 0.9565393
EN.ATM.CO2E.KT EN.ATM.NOXE.KT.CE 0.9550410
NE.IMP.GNFS.ZS NE.TRD.GNFS.ZS 0.9541148
SP.DYN.TO65.MA.ZS SP.DYN.TO65.FE.ZS 0.9529467
AG.LND.TOTL.K2 EN.ATM.GHGT.KT.CE 0.9520595
EN.ATM.NOXE.KT.CE SP.RUR.TOTL 0.9508327
AG.LND.TOTL.K2 SP.POP.TOTL.FE.IN 0.9497692
AG.LND.TOTL.K2 SP.POP.TOTL 0.9485983
AG.LND.TOTL.K2 SP.POP.TOTL.MA.IN 0.9469733
AG.LND.TOTL.K2 EN.ATM.METH.EG.KT.CE 0.9464401
NY.GDP.MKTP.KD.ZG NY.GDP.PCAP.KD.ZG 0.9458804
SP.POP.TOTL.FE.IN EN.ATM.GHGT.KT.CE 0.9456586
SP.DYN.CBRT.IN SP.POP.0014.TO.ZS 0.9454977
SP.POP.TOTL EN.ATM.GHGT.KT.CE 0.9443133
EN.ATM.CO2E.KT EN.ATM.METH.KT.CE 0.9441949
SP.POP.TOTL.MA.IN EN.ATM.GHGT.KT.CE 0.9426513
EN.ATM.CO2E.SF.KT EN.ATM.METH.KT.CE 0.9411552
AG.LND.TOTL.K2 SP.RUR.TOTL 0.9410677
SP.POP.0014.TO.ZS SP.POP.1564.TO.ZS -0.9393719
EN.ATM.CO2E.SF.KT EN.ATM.NOXE.KT.CE 0.9385839
EN.ATM.CO2E.LF.KT EN.ATM.NOXE.KT.CE 0.9355386
SP.RUR.TOTL SP.URB.TOTL 0.9341954
SP.DYN.LE00.IN SP.DYN.IMRT.IN -0.9322860
EN.ATM.CO2E.GF.KT SP.URB.TOTL 0.9306507
EN.ATM.CO2E.GF.KT EN.ATM.CO2E.SF.KT 0.9263083
EN.ATM.CO2E.SF.KT SP.POP.TOTL.FE.IN 0.9260594
EN.ATM.CO2E.GF.KT EN.ATM.METH.EG.KT.CE 0.9254724
EN.ATM.CO2E.SF.KT SP.POP.TOTL 0.9252594
EN.ATM.METH.EG.KT.CE SP.RUR.TOTL 0.9248062
EN.ATM.CO2E.SF.KT SP.POP.TOTL.MA.IN 0.9240043
AG.LND.TOTL.K2 SP.URB.TOTL 0.9230082
SP.DYN.IMRT.IN SP.DYN.TO65.FE.ZS -0.9164811
EN.ATM.CO2E.KT AG.LND.TOTL.K2 0.9164563
EN.ATM.CO2E.KT SP.POP.TOTL.FE.IN 0.9154624
EN.ATM.CO2E.GF.KT EN.ATM.NOXE.KT.CE 0.9140514
EN.ATM.CO2E.KT SP.POP.TOTL 0.9139805
EN.ATM.CO2E.KT SP.POP.TOTL.MA.IN 0.9120409
SH.DTH.MORT SP.RUR.TOTL 0.9119928
EN.ATM.CO2E.LF.KT EN.ATM.CO2E.SF.KT 0.9096632
EN.ATM.CO2E.GF.KT EN.ATM.METH.KT.CE 0.9091581
EN.ATM.CO2E.LF.KT EN.ATM.METH.KT.CE 0.9076110
EN.ATM.CO2E.LF.KT AG.LND.TOTL.K2 0.9068483
SP.DYN.CBRT.IN SP.POP.1564.TO.ZS -0.9048672
EN.ATM.CO2E.GF.KT NY.GDP.MKTP.CD 0.9047587
SP.RUR.TOTL EN.ATM.GHGT.KT.CE 0.8985380
EN.ATM.CO2E.LF.KT EN.ATM.METH.EG.KT.CE 0.8975581
EN.ATM.CO2E.LF.KT SP.URB.TOTL 0.8975480
EN.ATM.CO2E.SF.KT AG.LND.TOTL.K2 0.8970303
EN.ATM.CO2E.GF.KT AG.LND.TOTL.K2 0.8842807
AG.LND.TOTL.K2 SH.DTH.MORT 0.8828061
SP.DYN.CBRT.IN SP.DYN.TO65.FE.ZS -0.8822831
EN.ATM.CO2E.GF.KT SP.POP.TOTL.FE.IN 0.8775350
SP.DYN.CBRT.IN SP.DYN.LE00.IN -0.8770774
EN.ATM.CO2E.GF.KT SP.POP.TOTL 0.8757681
EN.ATM.CO2E.GF.KT SP.POP.TOTL.MA.IN 0.8732623
EN.ATM.CO2E.GF.KT NE.IMP.GNFS.CD 0.8678542
EN.ATM.CO2E.SF.KT SP.RUR.TOTL 0.8666477
EN.ATM.CO2E.GF.KT NE.EXP.GNFS.CD 0.8639633
EN.ATM.CO2E.LF.KT SP.POP.TOTL.FE.IN 0.8629516
EN.ATM.CO2E.LF.KT SP.POP.TOTL 0.8609962
EN.ATM.CO2E.LF.KT SP.POP.TOTL.MA.IN 0.8581721
SP.POP.0014.TO.ZS SP.POP.65UP.TO.ZS -0.8562991
SP.DYN.CBRT.IN SP.DYN.IMRT.IN 0.8562227
EN.ATM.CO2E.KT NY.GDP.MKTP.CD 0.8528268
SH.DTH.MORT SP.POP.TOTL.MA.IN 0.8524929
EN.ATM.CO2E.KT SP.RUR.TOTL 0.8524809
SH.DTH.MORT SP.POP.TOTL 0.8524546
SH.DTH.MORT SP.POP.TOTL.FE.IN 0.8519796
SP.DYN.IMRT.IN SP.DYN.TO65.MA.ZS -0.8518772
EN.ATM.METH.KT.CE SH.DTH.MORT 0.8515036
EN.ATM.NOXE.KT.CE SH.DTH.MORT 0.8498847
EN.ATM.CO2E.LF.KT NY.GDP.MKTP.CD 0.8275347
EN.ATM.CO2E.SF.KT NY.GDP.MKTP.CD 0.8118827
SP.POP.0014.TO.ZS SP.DYN.TO65.FE.ZS -0.8104836
NY.GDP.MKTP.CD EN.ATM.GHGT.KT.CE 0.8075622
EN.ATM.CO2E.LF.KT SP.RUR.TOTL 0.8060015
EN.ATM.CO2E.KT NE.IMP.GNFS.CD 0.8058971
EN.ATM.CO2E.GF.KT SP.RUR.TOTL 0.8056698
EN.ATM.CO2E.KT NE.EXP.GNFS.CD 0.8033274
SP.DYN.LE00.IN SP.POP.0014.TO.ZS -0.8011021
SP.DYN.CBRT.IN SP.DYN.TO65.MA.ZS -0.8000131
NY.GDP.MKTP.CD SP.URB.TOTL 0.7965653
SP.DYN.CBRT.IN SP.POP.65UP.TO.ZS -0.7844853
SP.POP.1564.TO.ZS SP.DYN.TO65.FE.ZS 0.7798112
EN.ATM.CO2E.LF.KT NE.IMP.GNFS.CD 0.7751604
EN.ATM.CO2E.LF.KT NE.EXP.GNFS.CD 0.7700196
SP.DYN.LE00.IN SP.POP.1564.TO.ZS 0.7695005
EN.ATM.CO2E.SF.KT NE.IMP.GNFS.CD 0.7678249
EN.ATM.CO2E.SF.KT NE.EXP.GNFS.CD 0.7675252
NE.IMP.GNFS.CD SP.URB.TOTL 0.7608635
NE.IMP.GNFS.CD EN.ATM.GHGT.KT.CE 0.7604903
SH.DTH.MORT EN.ATM.GHGT.KT.CE 0.7598616
NE.EXP.GNFS.CD SP.URB.TOTL 0.7594409
EN.ATM.METH.EG.KT.CE SH.DTH.MORT 0.7586380
NE.EXP.GNFS.CD EN.ATM.GHGT.KT.CE 0.7578417
SP.DYN.IMRT.IN SP.POP.0014.TO.ZS 0.7530084

 

Jak można było przewidzieć mamy tutaj do czynienia z wieloma korelacjami - z czego większość jest bardzo oczywistych, jak np. Population, female x Population, total czy CO2 emissions from gaseous fuel consumption (kt) x CO2 emissions (kt). Naszym zadaniem jest znalezienie najbardziej interesujących korelacji i jako takie najbardziej interesujące możemy wskazać:

  • Life expectancy at birth, total (years) oraz Survival to age 65, male (% of cohort) - wartość współczynnika korelacji Pearsona: 0,988
  • Methane emissions (kt of CO2 equivalent) oraz Number of under-five deaths - wartość współczynnika korelacji Pearsona: 0,852
  • GDP (current US$) oraz Population, total - wartość współczynnika korelacji Pearsona: 0,796

W przypadku pierwszej z korelacji mamy tutaj ciekawą zależność wraz ze wzrostem oczekiwanej długości życia danej osoby przy urodzeniu wzrasta wskaźnik dotyczacy przeżywalności 65 roku życia wśród mężczyzn. Świadczyć to może o nieustannej ewolucji i rozwoju medycyny - wraz z upływem lat metody służące do obliczania oczekiwanej długości życia zostają coraz bardziej dopracowane - dzięki rozwojowi medycznej części świata możliwości człowieka na dłuższe życie zwiększają się, co potwierdza znaleziona w danych korelacja.

Druga ze znalezionych korelacji jest również bardzo ciekawą przesłanką - wraz ze wzrostem emisji metanu wzrasta licba śmierci osób poniżej 5 roku życia. Korelacja ta ma najprawdopodobniej podłoże powiązane z populacją ludzi na świecie.

Jak możemy zauważyć na prezentowanym wykresie populacja ludzi na świecie notuje wzrost z roku na rok - wraz ze wzrostem populacji wzrasta też liczba śmierci osób poniżej 5 roku życia (naturalna konsekwencja wzrostu urodzeń). Dodatkowo wzrast ze wzrostem populacji wzrasta także emisja metanu - na świecie jest coraz więcej ludzi także również naturalnym wydaje się w tym przypadku wzrost emisji metanu, który jest głównym składnikiem gazu ziemnego.

Trzecia z korelacji jest ciekawym przypadkiem, który potwierdza wpływ wielkości populacji na wartość produktu krajowego brutto. Wraz ze wzrostem liczby ludności w danym kraju, można znaleźć coraz więcej rąk do pracy, co przekłada się na wzrost PKB mierzonego jako łączna wartość wszystkich dóbr i usług wytworzonym w danym kraju w ciągu roku.

Macierz wszystkich wyliczonych korelacji prezentuje się następująco:

Korelacje między cenami złota

W tej sekcji chcielibyśmy zbadać, które ze zmiennych są najbardziej skorelowane z wcześniej dołączaną do zbioru cechą zawierającą informację odnośnie ceny złota. Z racji, że wartości ceny złota dla każdego kraju są takie same, dlatego w tym przypadku zdecydowano się wybrać 4 kraje (USA, Polska, Indie, RPA), dla których zostaną policzone oraz wyszczególnione korelacje cech z cenami złota.

correlation country
NY.GSR.NFCY.CD 0.9453807 USA
NY.GSR.NFCY.CN 0.9453807 USA
NE.IMP.GNFS.CD2 0.9395134 Indie
NE.EXP.GNFS.CD2 0.9387731 Indie
NY.GDP.PCAP.CD2 0.9363200 Indie
NY.GDP.MKTP.CD2 0.9286484 Indie
NE.EXP.GNFS.CD1 0.9280784 Poland
NE.IMP.GNFS.CD1 0.9253574 Poland
NE.EXP.GNFS.CD3 0.9196490 RPA
NY.GDP.PCAP.CD1 0.9126322 Poland
NY.GDP.MKTP.CD1 0.9118309 Poland
NY.GDP.MKTP.CD3 0.9066883 RPA
EN.ATM.CO2E.KT 0.9063821 Indie
EN.ATM.GHGT.KT.CE 0.9048589 Indie
NE.EXP.GNFS.CD 0.8984441 USA
NE.IMP.GNFS.CD3 0.8923264 RPA
NE.TRD.GNFS.ZS 0.8848988 Poland
EN.URB.LCTY 0.8831704 RPA
EN.ATM.GHGT.KT.CE1 0.8826979 RPA
NY.GDP.PCAP.CD3 0.8791365 RPA
NE.IMP.GNFS.CD 0.8774351 USA
EN.ATM.CO2E.GF.KT 0.8681158 Poland
NY.GDP.MKTP.CD 0.8622578 USA
NY.GDP.PCAP.CD 0.8412912 USA

 

Jak możemy zauważyć w przypadku Polski, Indii i RPA najsilniejsze korelacje zostały znalezione w odniesieniu do takich cech związanymi z: wartością PKB, eksportem oraz importem dóbr i usług. W przypadku USA wyglądało to trochę odmiennie - największą korelację z ceną złota wykazały cechy związane z dochodami netto (Net primary income). W dalszej części postaramy się zweryfikować czy znalezione cechy rzeczywiście będą miały największy wpływ na predykcję cen złota przez stworzony regresor.

Badanie zmiany wybranych atrybutów w czasie

W tej sekcji zostaną zaprezentowane animacje wykresów prezentujące zmianę wybranych atrybutów w czasie.

S&P Composite

Wartości cen złota

Procentowy wzrost populacji a procentowy wzrost PKB

Na animacji być może słabo będzie widoczny moment spadku wartości PKB w roku 2020, ale jest to moment, na który warto zwrócić uwagę. Wartość spadła wtedy na około -4%, co ma swoje odzwierciedlenie w globalnym lock-down’ie spowodowanym przez pandemię wywołaną wirusem SARS-CoV-2. Warto także zwrócić uwagę na fakt, że oba te wskaźniki nie są ze sobą skorelowane jakby się mogło wydawać na początku - występują na animacji 4 następujące przypadki:

  • Wskaźnik PKB rośnie, wskaźnik populacji maleje
  • Wskaźnik PKB rośnie, wskaźnik populacji rośnie
  • Wskaźnik PKB maleje, wskaźnik populacji maleje
  • Wskaźnik PKB maleje, wskaźnik populacji rośnie

Problem predykcji cen złota

W tej sekcji zostanie zaprezentowana próba stworzenia regresora przewidująca ceny złota. Wybrano do predykcji ceny złota ze względu na znacznie większą ilość danych posiadanych (od 1970 do 2020 roku) aniżeli dla Bitcoina, gdzie dane posiadamy dopiero od roku 2009.

Przewidywanie cen złota z wykorzystaniem zbioru World_Development_Indicators

W pierwszej próbie wykorzystamy wcześniej utworzony i analizowany zbiór danych bazujący głównie na danych ze zbioru World_Development_Indicators. Głównym problemem w tym podejściu jest fakt, że wartości cen złota z dokładności do jednego dnia zostały sprowadzone do jednej wartości na rok. Drugim problemem jest fakt, że dla różnych krajów, ale dla określonego roku ceny złota powtarzają się - globalna cena złota taka sama dla każdego z krajów. W obliczu opisanych problemów możemy zastosować dwa podejścia:

  • Przeprowadzenie predykcji dla danych dla wybranego kraju
  • Przeprowadzenie predykcji dla danych dostępnych dla wszystkich krajów

W pierwszym podejściu do procesu predykcji zachowalibyśmy jedynie 52 obserwacje - jest to bardzo mało dla procesu uczenia maszynowego przez co końcowa predykcja mogłaby wypaść bardzo słabo. W związku na przedstawioną wadę pierwszego podejścia zastosujemy do procesu uczenia cały dostępny zbiór danych - z drobnymi wyjątkami, o których mowa w następnej sekcji.

Przygotowanie zbioru danych do procesu uczenia

W utworzonym przez nas zbiorze danych zostaną poczynione małe modyfikacje, mianowicie - usunięte zostaną kolumny zawierające ceny Bitcoina, S&P Composite oraz Real Price, z tego względu, że tak naprawdę klasyfikator byłby w stanie przewidzieć ceny złota tylko na podstawie tych zmiennych, z tego samego powodu, o którym napisano w poprzedniej sekcji - wartości z całego roku zostały zebrane do jednej wartości i powtarzają się dla danego roku niezależnie od kraju.

train_df <- clean_df %>%
  select(-c(BC_PRICE, COMPOSITE, PRICE, USD_AM))

train_df <- train_df[complete.cases(train_df),]

Podział zbioru danych dla procesu uczenia

set.seed(42)
inTraining <- createDataPartition(y = train_df$USD_PM, p = .75, list = FALSE)

training <- train_df[inTraining,]
testing <- train_df[-inTraining,] %>%
  filter(c_code %in% unique(training$c_code))

Pierwsze 75% zbioru trafiło do zbioru treningowego, natomiast pozostała reszta znalazła się w zbiorze testowym. Do podziału zbioru skorzystano ze standardowej metody createDataPartition.

Proces uczenia z wykorzystaniem algorytmu Random Forest

Jako algorytm predykcyjny zostanie wykorzystanu algorytm Random Forest, który jest absolutnym algorytmem podstawowym, który warto zaaplikować jako pierwszy do swoich danych jako rozwiązanie podstawowe (baseline). W trenowaniu zostanie wykorzystana walidacja krzyżowa z podziałem na 5 części.

rfGrid <- expand.grid(mtry = 20:30)
gridCtrl <- trainControl(method = "cv", number = 5)

set.seed(42)
rfFitTune <- train(USD_PM ~ .,
                   data = training,
                   method = "rf",
                   trControl = gridCtrl,
                   tuneGrid = rfGrid,
                   ntree = 50)

rfFitTune
## Random Forest 
## 
## 5242 samples
##   54 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 4193, 4195, 4194, 4194, 4192 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   20    120.0531  0.9402284  80.87823
##   21    120.3805  0.9399695  80.77626
##   22    116.8532  0.9430249  78.38126
##   23    116.4319  0.9438273  78.23815
##   24    114.5341  0.9456541  77.20251
##   25    111.4655  0.9477998  75.23970
##   26    109.3990  0.9499193  73.87188
##   27    109.0996  0.9502753  73.84088
##   28    104.4382  0.9539389  70.65630
##   29    103.1778  0.9556620  70.22044
##   30    101.7334  0.9562226  68.54536
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 30.

Predykcja

Predykcja modelu zostanie przebadana pod kątem trzech podstawowych metryk: RMSE, MAE i R^2 szczególnie zwracając uwagę na tę pierwszą metrykę, z tego względu, że RMSE w stosunku do MAE zdecydowanie bardziej zwraca uwagę na wartości odstające (outliery) i przypisuje im większą wartość aniżeli MAE. Z kolei RMSE w stosunku do R^2 zdecydowanie lepiej radzi sobie z wykryciem problemu przeuczenia modelu, czego nie można powiedzieć o metryce R^2.

rfTuneClasses <- predict(rfFitTune, newdata = testing)
RMSE(rfTuneClasses, testing$USD_PM)
## [1] 94.69507
MAE(rfTuneClasses, testing$USD_PM)
## [1] 64.07353
rsq <- function(x, y) {
  cor(x, y) ^ 2
}

rsq(rfTuneClasses, testing$USD_PM)
## [1] 0.9611175

Otrzymane wyniki metryk prezentują się następująco:

  • RMSE = 94,695
  • MAE = 64,07
  • R^2 = 0,961

Jak możemy zauważyć z otrzymanej wartości R^2 możemy stwierdzić, że 96% zmienności może zostać wyjaśnione przez nasz model. Na podstawie wartości RMSE możemy stwierdzić, że średnia wartość ceny złota między obserwowanymi wartościami danych a przewidywanymi wartościami danych wynosi lekko ponad 94 dolary.

Przewidywanie cen złota z pominięciem zbioru World_Development_Indicators

W celu poprawy otrzymanego w poprzedniej sekcji wyniku postanowiono sprawdzić jak nowo stworzony model poradzi sobie z problemem predykcji cen złota. W tym przypadku nie skorzystamy ze zbioru World_Development_Indicators, a naszą predykcję przeprowadzimy w oparciu o pozostałe zbiory danych, gdzie posiadamy dane z szczegółowością do jednego dnia aniżeli do jednego roku. Sprawdzimy w ten sposób czy zapewnienie większej szczegółowości danych wpłynie na lepsze wyniki modelu aniżeli dla przypadku, gdy cech dla danej obserwacji jest znacznie więcej i szczegółowość czasu jest mniejsza.

Przygotowany zbiór danych po połączeniu i uzupełnieniu wartości pustych prezentuje się następująco

##         Date USD_PM BC_PRICE COMPOSITE  CPI RATE BC_DIFF BC_HRATE
## 1 1968-01-02   37.7        0     95.04 34.1 5.53       1        0
## 2 1968-01-03   37.7        0     95.04 34.1 5.53       1        0
## 3 1968-01-04   37.7        0     95.04 34.1 5.53       1        0
## 4 1968-01-05   37.7        0     95.04 34.1 5.53       1        0
## 5 1968-01-08   37.7        0     95.04 34.1 5.53       1        0
## 6 1968-01-09   37.7        0     95.04 34.1 5.53       1        0

Podział zbioru danych dla procesu uczenia

set.seed(42)
inTraining_sec <- createDataPartition(y = train_sec_df$USD_PM, p = .75, list = FALSE)

training_sec <- train_sec_df[inTraining_sec,]
testing_sec <- train_sec_df[-inTraining_sec,]

Pierwsze 75% zbioru trafiło do zbioru treningowego, natomiast pozostała reszta znalazła się w zbiorze testowym. Do podziału zbioru skorzystano ze standardowej metody createDataPartition.

Proces uczenia z wykorzystaniem algorytmu Random Forest

rfGrid_sec <- expand.grid(mtry = 20:30)
gridCtrl_sec <- trainControl(method = "cv", number = 5)

set.seed(42)
rfFitTune_sec <- train(USD_PM ~ .,
                   data = training_sec,
                   method = "rf",
                   trControl = gridCtrl_sec,
                   tuneGrid = rfGrid_sec,
                   ntree = 50)

rfFitTune_sec
## Random Forest 
## 
## 10190 samples
##     7 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8153, 8153, 8151, 8152, 8151 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   20    9.825482  0.9995966  5.092764
##   21    9.824705  0.9995961  5.130806
##   22    9.916169  0.9995882  5.122228
##   23    9.778131  0.9995995  5.099669
##   24    9.873135  0.9995913  5.095654
##   25    9.855161  0.9995933  5.120413
##   26    9.882656  0.9995909  5.135560
##   27    9.900777  0.9995900  5.142456
##   28    9.970601  0.9995842  5.155781
##   29    9.812591  0.9995980  5.124667
##   30    9.826999  0.9995960  5.130117
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 23.

Predykcja

rfTuneClasses_sec <- predict(rfFitTune_sec, newdata = testing_sec)
RMSE(rfTuneClasses_sec, testing_sec$USD_PM)
## [1] 8.567937
MAE(rfTuneClasses_sec, testing_sec$USD_PM)
## [1] 4.495058
rsq(rfTuneClasses_sec, testing_sec$USD_PM)
## [1] 0.9996967

Otrzymane wyniki metryk prezentują się następująco:

  • RMSE = 8,57
  • MAE = 4,49
  • R^2 = 0,999

Na podstawie otrzymanych wyników możemy stwierdzić, że drugi stworzony przez nas model prezentuje się lepiej aniżeli pierwszy model, co potwierdzają wyniki poszczególnych metryk. Na podstawie RMSE możemy stwierdzić, że średnia wartość ceny złota między obserwowanymi wartościami danych a przewidywanymi wartościami danych wynosi lekko ponad 8,5 dolara. Około 99% zmienności może zostać wyjaśnione przez nasz zbudowany model - metryka R^2.

Analiza ważności atrybutów

Analiza ważności atrybutów zostanie przeprowadzona dla najlepszego znalezionego modelu.

varImp(rfFitTune_sec)
## rf variable importance
## 
##             Overall
## Date      1.000e+02
## CPI       4.900e+00
## RATE      4.325e+00
## COMPOSITE 2.989e-01
## BC_DIFF   3.819e-02
## BC_PRICE  1.428e-03
## BC_HRATE  0.000e+00

Jak możemy zauważyć na podstawie otrzymanych wyników widzimy, że największe znaczenie w procesie predykcji ma data odnotowania wpisu wartości złota, co ciekawe atrybut BC_HRATE nie jest w ogóle brany pod uwagę w procesie predykcji. Warto zwrócić uwagę na fakt, że ważność między datą a drugim atrybutem pod tym względem (CPI) jest ogromna. Skąd może wynikać taka różnica? Mianowicie głównym problemem w tym przypadku jest sposób losowania przykładów do zbioru treningowego i testowego. Mamy tutaj do czynienia z danymi w kontekście czasu - wykorzystana metoda createDataPartition nie zwraca uwagi na kolumnę zawierającą aspekt czasowy tylko dobiera przykłady w ten sposób, aby rozkłady cen złota w zbiorze treningowym i testowym były jak najbardziej podobne, w efekcie czego nastąpiło przemieszanie obserwacji między zbiorami. Przykładowo w zbiorze treningowym znalazły się wpisy z dni 10.10.2020 oraz 12.10.2020, a z zbiorze testowym wpis z dnia 11.10.2020 w efekcie czego model miał bardzo ułatwione zadanie, gdyż znał wartości ceny złota w bliskim okresie od daty obserwacji i mógł sprawnie podać wartość ceny złota na dzień 11.10.2020, która niewiele różniła się od rzeczywistej ceny, gdyż wartości cen złota nie notują bardzo dużych spadków w pojedynczy dzień.

W jaki sposób moglibyśmy sobie poradzić z problemem predykcji wyłącznie za pomocą atrybutu daty? Można by było zastosować inną techniką tworzenia zbioru testowego i zbioru treningowego - mianowicie pierwsze 75% przypadków trafiłoby do zbioru treningowego, a pozostała część do testowego. W ten sposób uniknęlibyśmy zjawiska przemieszania przykładów między zbiorami. Niestety pojawia się jednak kolejny problem,

Jak spojrzymy na wykres cen złota możemy zauważyć, że od około roku 2005 zaczął się szybki wzrost cen złota, który trwał aż do roku około 2013. Przy założeniu, że zbiór treningowy zawierałby 75% całego zbioru to jego ostatnia obserwacja wystąpiłaby w roku 2008. Oznacza to, że w zbiorze testowym znaczna większość obserwacji dotyczyłaby trendów szybkiego wzrostu oraz szybkiego spadku cen, podczas gdy w zbiorze treningowym bardzo rzadko występowały sytuacje szybkich spadków lub wzrostów. Stworzony na takich danych model mógłby sobie po prostu nie poradzić z charakterystyką danych znajdujących się w zbiorze testowym, gdyż uczony był na zupełnie innej charakterystyce.

W celach porównawczych został stworzony właśnie taki model danych i tak jak się spodziewano otrzymane wyniki były bardzo słabe:

  • RMSE = 530,12
  • MAE = 461,46
  • R^2 = 0,12

Analiza aspektu czasu z wykorzystaniem metody ARIMA

Jako ostatni element naszej analizy sprawdzimy jak poradzi sobie z naszymi danymi dedykowana metoda do analizy szeregów czasowych - ARIMA

Zbudowanie zbioru treningowego oraz testowego

train_range <- 1:(0.75 * nrow(train_sec_df))

train_arima_ts <- ts(train_sec_df[train_range,])
test_arima_df <- train_sec_df[-train_range,]

Pierwsze 75% obserwacji trafiło do zbioru treningowego, pozostała częśc trafiła do zbioru testowego.

Proces uczenia z wykorzystaniem algorytmu ARIMA

set.seed(42)
arimaFitTune <- auto.arima(train_arima_ts[, 2])

arimaFitTune
## Series: train_arima_ts[, 2] 
## ARIMA(4,1,2) with drift 
## 
## Coefficients:
##          ar1      ar2     ar3      ar4      ma1     ma2   drift
##       1.1318  -0.6041  0.0517  -0.0902  -1.1711  0.6365  0.0863
## s.e.  0.0693   0.0594  0.0163   0.0107   0.0695  0.0601  0.0467
## 
## sigma^2 estimated as 26.74:  log likelihood=-31190.07
## AIC=62396.13   AICc=62396.14   BIC=62453.96

Predykcja

Tak prezentują się predykowane wartości przez algorytm ARIMA w porównaniu do wartości rzeczywistych cen złota.

accuracy(forecast(arimaFitTune, h = nrow(test_arima_df)), test_arima_df[, 2]) %>%
  kable %>%
  kable_styling(c("striped", "hover"))
ME RMSE MAE MPE MAPE MASE ACF1
Training set 0.0000187 5.169424 2.692902 -0.0306804 0.7911418 1.001665 0.0005225
Test set 285.9467123 373.968125 296.588859 18.5322466 19.8383704 110.320616 NA

 

Jak możemy zauważyć wyniki otrzymane z wykorzystaniem algorytmu ARIMA dla przedstawionej generacji zbioru treningowego i testowego są lepsze aniżeli przy wykorzystaniu standardowego algorytmu Random Forest. Wartość metryki RMSE spadła o prawie 200 jednostek, wartość metryki MAE spadła o około 170 jednostek. Prowadzi nas to do wniosku, że warto w danych charakteryzujących się przebiegiem czasowym skorzystać z dedykowanych metod do analizy, ponieważ jest duża szansa, że nasze wyniki predykcji ulegną poprawie tak samo jak to się stało w naszym eksperymencie.